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Abstract 

A Lennard-Jones model of a binary dense liquid (A,B) with a symmetrical miscibility gap is 
investigated by means of computer simulation methods. Semigrand-canonical Monte Carlo sim- 
ulations yield the phase diagram in the T—x plane (T: temperature, x: concentration of A or B 
particles) as well as equilibrated configurations at coexistence. Then Molecular Dynamics simu- 
lations use these configurations to determine static properties (isothermal compressibility kt and 
concentration susceptibility x) as wen as the shear and the bulk viscosity rj B and r/B, respectively. 
The latter quantities are calculated along a path approaching the coexistence line from high tem- 
peratures in the one-phase region and ending at a state at the coexistence line about 15% below 
the critical point. We find that kt and x increase significantly near the coexistence line reflecting 
the vicinity of the critical point. Whereas % exhibits a weak temperature dependence, t]b increases 
significantly near the coexistence curve. 
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I. INTRODUCTION 



The bulk viscosity (in the following denoted by t/b) describes the response of a fluid to 
a compression or expansion. Compared to other transport coefficients such as the shear 
viscosity or the self diffusion constant, it is the least studied transport coefficient. This is 
surprising since t/b is for instance a central quantity in the description of the damping of 
longitudinal sound. It is also an important quantity to probe slow dynamic processes such as 
the critical slowing down near the critical point of a liquid-gas transition or the liquid-liquid 
unmixing transition in a binary fluid. We will briefly discuss these issues below. 

A microscopic expression for rj^ is given by a Green-Kubo formula (Boon and Yip, 1980), 

V r°° 

r lB = — (J aa (t)J aa (0)) , (I) 
k B l Jo 

with a denoting Cartesian components (a G {x, y, z}). V, T and ks are volume, temperature 
and Boltzmann's constant, respectively. In the microcanonical ensemble, J aa is equal to the 
difference between the pressure at time t and that at t = 0, J aa (t) = p(t) — p(0), where p(t) 
is equal to the diagonal elements of the pressure tensor a defined as follows: 
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Herein miVi a and ri a are respectively the a'th component of momentum and position of 
particle i, and Fi a is the a'th component of the force acting on particle i. Note that in order 
to calculate the shear viscosity ?7 S one has to use the non-diagonal elements of the pressure 
tensor in the Green-Kubo integral (a ^ (3) (Boon and Yip, 1980): 

V f°° 

^ = TT^r / dt (^(<)^(0)) • (3) 

Kgl JO 

Eqs. (JTJ), (J2J), and (JHJ) can be used to calculate t/b and r] s from equilibrium fluctuations in 
a Molecular Dynamics (MD) computer simulation. Indeed, in one of the pioneering MD 
studies of a Lennard- Jones liquid near its triple point by Levesque et al. (Levesque et 
al., 1973) the viscosities were determined by Green-Kubo formulas. Note that a recently 
proposed "new" formula by Okumura and Yonezawa (Okumura and Yonezawa, 2002) just 
expresses the pressure fluctuations in Eq. (JT} in terms of the pair correlation function and 
the interatomic potentials. 

Alternative methods to determine t]b are based on Non-Equilibrium Molecular Dynamics 
(NEMD) simulations. Heyes (Heyes, 1984; Heyes, 1986) proposed a NEMD scheme where 
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the volume of the system is changed from V to V + AV at t = which leads to a change of 
the pressure. Then one follows the relaxation of the pressure to its equilibrium value p(oo) 
and measures p(t) — p(oo) where p(t) is the instantaneous pressure at time t. r/ B is then 
given by 



Certainly, Eq. is only valid if AV is small enough to allow the application of linear 
response theory. 

Another NEMD approach was proposed by Hoover et al. (Hoover et al., 1980). The latter 
authors impose a frequency-dependent small perturbation by changing the volume of the 
system by a periodic compression and expansion with a frequency u. As a result ^(oj) is 
obtained for several values of uj and then an extrapolation to zero frequency may be possible. 
The method of Hoover et al. might be especially useful for liquid states where t/b exhibits a 
long-time tail. However, some knowledge of the frequency dependence of r] B (u) is required 
to extrapolate it accurately from finite frequencies to zero. 

The aforementioned methods have been mainly used in feasibility studies where the bulk 
viscosity was determined, e.g., for a Lennard- Jones fluid at a single state near its triple point 
(see Refs. (Levesque et al., 1973; Hoover et al., 1980; Heyes, 1984; Okumura and Yonezawa, 
2002)). Only in a small number of simulations, t]b has been investigated systematically. 
One of these rare studies is the MD simulation of symmetrical Lennard- Jones mixtures 
by Vogelsang and Hoheisel (Vogelsang and Hoheisel, 1988; Hoheisel, 1993) who considered 
systems of 256 particles at moderate densities (i.e. far from the triple point). In this work 
t)b as well as rj s were calculated by means of the Green-Kubo formulas, Eqs. (0) and (J3J. An 
interesting result of this study was that the ratio t/b/^s is (much) larger than one if the fluid 
mixture has a (strongly) associating character or a (strongly) demixing character. In both 
of the latter cases the bulk viscosity increases quickly whereas the shear viscosity remains 
essentially constant. As a consequence it is expected that r] B shows a strong increase near 
the coexistence line of a fluid-fluid unmixing transition. 

In contrast to the small number of simulations, there are many theoretical investigations of 
the bulk viscosity in the context of the dynamics near the liquid-gas critical point (Kawasaki, 
1976; Kadanoff and Swift, 1968; Swift, 1968; Hohenberg and Halperin, 1977; Folk and Moser, 
1995; Onuki, 1997; Onuki, 2002). These works predict that the bulk viscosity exhibits a 
strong divergence near the critical point of a gas-liquid liquid transition. In contrast to that, 
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the shear viscosity is expected to show a very weak divergence (logarithmic divergence) at 
the critical point (if there is at all a divergence in this quantity). The latter predictions have 
been confirmed experimentally. An example is 3 He in the vicinity of the critical temperature 
T c : At T/T c — 1 = 10~ 4 on the critical isochore, rfe is about 50 Poise whereas ?7 S is equal to 
17 x KT 6 Poise (Kogan and Meyer, 1998; Onuki, 2002). 

In the present work we consider a simple model of a dense liquid mixture near and 
at a liquid-liquid unmixing transition and, apart from static susceptibilities, we calculate 
the shear and the bulk viscosity. Although we are not able to determine these quantities 
very close to the critical point, we find a behavior which agrees qualitatively with the 
aforementioned theoretical predictions for the critical dynamics: t]-q shows a stronger increase 
than ?7 S when approaching a state on the coexistence line about 15% below the critical point 
and, furthermore, at the latter point, tjb is significantly larger than ?7 S , i.e. r/e/^s ~ 3.3. 

The rest of the paper is organized as follows: In the next section we briefly comment 
on the details of the simulation as well as the Lennard- Jones model and its phase diagram. 
The static properties and the transport coefficients (shear and bulk viscosity) as obtained 
from the simulation are then presented in Sec. 3. Finally we summarize the results in Sec. 4. 

II. MODEL AND PHASE DIAGRAM 

The model that we consider in this work is a binary Lennard- Jones mixture. Thus, the 
interaction potential between a particle of type a and a particle of type j3 (a, f3 G {A, B}) 
is given by 



r being the distance between the two particles. For the Lennard-Jones parameters e a p and 
<j a /3 we choose <taa = obb = °~ab = cr, caa = £bb = e and cab = Se. Lengths, energies, 
and temperatures are measured respectively in units of a = 1, e = 1, and e/ks = 1. 
In the Molecular Dynamics (MD) part equal masses are chosen for A and B particles, 
i.e. mA = m-Q = 1. The potential is truncated and shifted at r = 2.5o\ 

The model mixture that we have defined so far is obviously completely symmetrical. 
Whether it has the tendency towards association or demixing is controlled by the parameter 
5. We use 5 = 0.5 which implies the possibility of a fluid-fluid unmixing transition. Since 
we are interested in the dense liquid state we have chosen a density pa 3 = 1, which provides 
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the absence of crystallization in the temperature range of interest, T > 1.0. Note that for 
densities < p < 0.7 the phase behavior of symmetrical LJ mixtures have been extensively 
investigated by Wilding (Wilding, 1997). 

The simulations were done as follows: We started with a random mixture with an equal 
number of A and B particles. By using standard Monte Carlo (MC) in the canonical ensemble 
with trial displacements of particles in the range [— a/20, +<r/20], we equilibrated the system 
for 10 5 Monte Carlo steps (MCS) per particle. Then, we switched on a MC simulation in 
the semigrand-canonical ensemble, i.e., at the end of each displacement sweep an identity 
switch of iV/10 randomly chosen particles was attempted, A — > B or B — > A (N being the 
total number of particles). Note that in the Metropolis criterion of the semigrand-canonical 
moves, the chemical potential energy difference ±A/z = — /i B (p a - chemical potential of 
species a G {A, B}) has to be taken into account in addition to the energy change in the 
Boltzmann factor. In order to localize the coexistence curve of the liquid-liquid unmixing 
transition in the present case, one has to just set A/i = which is simply due to the symmetry 
of our model. In order to determine the phase diagram we have performed five independent 
runs with a length of 400000 MCS per particle where we started the averaging after 100000 
MCS in each run (for more details of this calculation see Ref. (Das et al., 2003)). 

Fig.^shows the phase diagram in the T-xb plane for the system sizes N = 400, 800, 1600, 
and 3200 (xb = N^/N is the concentration of B particles). Due to the symmetry of the 
model we know a priori that the critical point is located at xb = 0.5. As we can infer from 
Fig. ^ the finite size effects near the critical point are small for N > 400, and for iV > 1600 
the data agree within the statistical errors. We have estimated the critical temperature T c 
from power law fits according to the three-dimensional Ising universality class (Binder and 
Ciccotti, 1996; Binder and Heermann, 2002), 

f(x B ) = 0.5 ± x c B ocx = B (1 - T/T c f , (3 « 0.325 (6) 

where B is a critical amplitude which is used, as well as T c , as a fitting parameter. From the 
fits with Eq. (jUJ) we obtain T c rs 1.638 ± 0.005 for N > 1600. For a more accurate estimate 
of T c , we would have to perform a finite scaling analysis (Binder and Ciccotti, 1996; Binder 
and Heermann, 2002). 

Apart from the phase diagram the MC in the semigrand-canonical ensemble yields also 
equilibrated configurations exactly along the coexistence line. We used them as starting con- 
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figurations for Molecular Dynamics (MD) simulations to determine the static quantities and 
the transport coefficients that are presented in the next section. In the MD, the equations 
of motion were integrated by means of the velocity Verlet algorithm (Binder and Ciccotti, 
1996) with a time step 5t = 0.01 [in units of the time t = (mcx 2 /( 48e )) 1/2 ]- 

Starting point for the MD were the configurations with 1600 particles at T = 1.4 that 
correspond to the concentration xb = 0.10375 at the coexistence line. Configurations in 
the one-phase region at the latter value of xb were obtained by heating up the system and 
equilibrating it for 10 5 time steps at constant temperature with the use of an Andersen 
thermostat (Das et al., 2003). Then, microcanonical runs were added from which we com- 
puted the static and dynamic quantities that are shown in the next section. The path along 
which we determined the latter quantities is indicated in Fig. by crosses: Apart from the 
coexistence state at T = 1.4, which is about 15% below the critical point with respect to 
temperature, the temperatures T = 1.6, 1.7, 3.0, and 6.0 were analyzed (note that also other 
paths around the coexistence line are studied in Ref. (Das et al., 2003)). 

One may wonder why we have not studied states that are much closer to the critical 
point. But due to the diverging correlation length that is accompanied by the approach of 
the critical point, we would have to consider systems that contain much more than 1600 
particles as in our work. Furthermore, the critical slowing down would require very long 
runs to equilibrate the system and to determine the transport coefficients with reasonable 
error bars. The latter point is especially a severe problem for transport coefficients such as 
the shear or the bulk viscosity. These are collective quantities and require many independent 
runs and/or a long time averaging since they are not subject to a self-averaging as one- 
particle quantities such as the self-diffusion constant. However, compared to many previous 
works, our choice of N is relatively large. Even the very recent computation of the bulk 
viscosity by Okumura and Yonezawa (Okumura and Yonezawa, 2002) was only done for a 
small system of 256 particles. 

III. RESULTS 

In this section we present the results for the static and dynamic properties of the sym- 
metrical LJ system along the path in the phase diagram which is indicated in Fig. ^ As 
described in the previous section, we have generated first five independent configurations 
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at each temperature. All these configurations were used as initial configurations for mi- 
crocanonical MD runs over 4.8 million time steps (the time step was 0.01 to, see previous 
section). Thus, at each temperature 24 million time steps were done to determine the quan- 
tities of interest. As we shall see in the following, this effort was large enough to get a 
reasonable estimate of bulk and shear viscosities. 



A. Static Properties 

As we see in Fig. Uthe states at T = 1.4 on the coexistence line are about 15% below 
the critical point with respect to temperature. Although these points are not very close 
to the critical point one may expect that the approach of the critical point is reflected in 
thermodynamic quantities such as the isothermal compressibility and the concentration 
susceptibility x (defined below). 

k,t can be calculated from the static number-number density structure factor S nn (q) in 
the limit of wavenumber q — > (Hansen and McDonald, 1986), 

«r = -7^1imS nn (g) (7) 

with p being the total density of the system (in our case p as well as the Boltzmann constant 
ks are equal to one). The structure factor S nn (q) for a binary AB mixture is defined 
by (Hansen and McDonald, 1986) 

S im (q) = S AA (q) + 2S AB (q) + S BB (q) (8) 

where S a p(q) (a, (3 G [A, B]) are the partial structure factors, 



M?) = ^E(^p[iq-(r?-rf)j) (9) 

with f a p = 0.5 for a ^ (3 and f a/ 3 = 1.0 for a = (3. In Eq. (JHJ) the indices i,j run over 
the number of particles of species a and (3, respectively, and rf is the position of the z'th 
particle of species a. 

Fig. |21 shows S nn (q) for T = 1.4, 1.7, 3.0, and 6.0. For wavenumbers q that correspond 
to distances smaller or equal the typical nearest neighbor distance, say q > 5, the typical 
behavior of this quantity for simple dense liquid can be identified: Upon decreasing the 
temperature the amplitude especially of the first peak increases and the peaks become 



narrower. The small values of S nn (q) for q — > reflect the fact that the considered dense 
liquid state is hardly compressible. It might be surprising that even at coexistence S nn does 
not show any tendency to increase significantly for q — > 0. The amplitude of S nn (q) at 
small q appears to be even a monotonic decreasing function with decreasing temperature. 
However, the relevant thermodynamic quantity in our context is Kt, that we have extracted 
from S nn (q) by extrapolating this function to q = 0. As we see in the inset of Fig. |2l 
kt increases significantly with decreasing temperature which shows that for states around 
T = 1.4, long-ranged static correlations, i.e. the presence of the critical point, still affect 
the behavior of Kt- 

The "vicinity" of the critical point is more pronounced in the structure factor of the 
concentration densities, S cc (q), than in S nn (q). S cc (q) can be also expressed by a linear 
combination of the partial structure factors (Hansen and McDonald, 1986), i.e. 

S cc {q) = x B S AA (q) - 2x A x B S AB (q) + x A S BB (q) . (10) 

In the limit q — > the structure factor S cc (q) is related to the static concentration suscepti- 
bility x by 

X = — ™ h m S cc (q) . (11) 

Note that we have determined \ directly via a fluctuation relation by semigrand-canonical 
MC runs (see Ref. (Das et al., 2003)). So, it was not necessary to extrapolate S cc (q) to 
q = 0. As we see in Fig. El this would be a difficult task because, in contrast to S nn (q), 
S cc (q) steeply increases for q — > 0. As we can infer from the inset of Fig. El X is about a 
factor of 2 larger at the coexistence state at T = 1.4 than at T = 1.8. It is remarkable that 
S cc (q) exhibits almost no temperature dependence for q > 5 in the broad temperature range 
1.4 < T < 6.0. 

B. Bulk viscosity and shear viscosity 

For the computation of the bulk and shear viscosities we have used the Green-Kubo 
(GK) formulas, Eqs. (0) and (JHJ)- The alternative methods that are based on NEMD require 
essentially the same computational effort. Furthermore, in the Heyes method, see Eq. (0), 
one has to choose the perturbation AV small enough to ensure that this perturbation justi- 
fies the application of linear response theory. Thus, one has to study the dependence of the 
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measured bulk viscosity t/b on AV^ (of course, in the linear response regime the apparent tjb 
is independent of AV). The Hoover method has in addition the drawback that one has to 
extrapolate the frequency-dependent viscosity to zero frequency. However, a comparative 
study of the different NEMD and GK methods to measure transport coefficients in a simula- 
tion is an interesting future project since the NEMD methods may give additional physical 
insight into the microscopic mechanism of different transport processes. 

Fig. |U shows r) s (t) and 773(2) for four temperatures. These quantities are defined by 
Eqs. (JTJ) and (J3J) where one has to replace oo in the integral by 2. We see that r] s (t) and 773(2) 
approach indeed plateaus at long times the values of which correspond to the hydrodynamic 
shear and bulk viscosities, respectively. At low temperatures, there is a qualitative difference 
in 7/b(2) as compared to i] s (t): E.g. at T = 1.4, 77 s (2) is essentially constant for 2 > 10. In 
contrast to that, 773(2) exhibits a second strong increase and it reaches the plateau value 
for 2 > 300. This is due to a long-time tail in the autocorrelation function of the pressure 
fluctuations. Note that the decrease of 773(2) for 2 > 500 is due to the fact that the statistics 
is much worse at long times. 

77B and 77s are plotted in Fig.GOa as a function of inverse temperature. Whereas r] s exhibits 
only a very weak temperature dependence, 773 increases significantly in the vicinity of the 
coexistence state at T = 1.4. As we can see in Fig. ISJd the ratio 773/778 is in the whole 
considered temperature range 6.0 > T > 1.4 above one, and it reaches a value of about 
3.3 at T = 1.4. One expects such a behavior from theories of the critical dynamics of the 
liquid-gas transition (Onuki, 2002). According to these theories the long-ranged critical 
fluctuations cause a slowing down of the system's response to a compression or expansion 
(described by 773). On the other hand, the response to the shearing of the system is hardly 
affected by the critical fluctuations (and thus 77s). In our case, at a state about 15% below 
the critical point, there is already a significant increase of static correlations which makes 
the behavior of 77b/t7 s plausible. 

Since the data presented in this paper are taken at an off-critical concentration, one could 
also attempt to interpret them in terms of a singular behavior at the "spinodal temperature" 
T s (limit of metastability) (Binder, 1987). According to the mean field theory of symmetric 
binary mixtures, one should expect that the static concentration susceptibility x f° r x b < 
x B nt = 0.5 behaves as 

X (T,x B )oc[T-T s (a:Br 1 (12) 
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where near the critical temperature the spinodal temperature T s (xb) is the inverse function 
of the concentration x B (T) along the spinodal curve, given by the equation x B (T) — x B nt = 
(:r B oex -x B rit )/^3 (Binder, 1987). Further away from T c , a simple expression for T s (xb) exists 
for the lattice (Ising) model of symmetric binary mixtures, namely 



Here we have emphasized by this notation that the mean field estimate T c of the critical 
temperature for systems with short range forces normally exceeds the actual critical tem- 
perature distinctly (also Eq. ()12j) does then not hold for x-q near ig 1 ' and T near the actual 
critical temperature, since x(T, i| nt ) oc (T — T c )~ 7 , where the actual susceptibility exponent 
7 « 1.24 (Binder and Ciccotti, 1996; Binder and Heermann, 2002)). 

Although we do not really expect that Eq. ()12)) and a related mean-field divergence for 
the bulk viscosity i]b is a good approximation for our Lennard- Jones system, we present a 
plot of vs. T and 773 1 vs. T in Fig El Mean field theory would predict that the data 
fall on straight lines and both straight lines should hit the abscissa in the same point which 
then is the estimate of T s (x-b)- Indeed the data points close to the coexistence curve are 
compatible with such analysis, with T s (xb) ~ 1. Of course, one should not put too much 
weight on this analysis, since the temperature range over which we need to extrapolate is 
larger than the temperature range where actual data are fitted. Also the estimate from 
Eq. (fT3*j) would be much lesser, namely T s (xb) — 0.59, if the distinction between the actual 
T c and T C MF is ignored. We caution the reader that anyway the concept of a spinodal is of 
doubtful validity outside of mean field theory (Binder, 1987), although in the experimental 
literature on binary mixtures (both in metallic alloys and in polymer blends, for instance) 
it is widely used. 

IV. SUMMARY 

We have used computer simulations to investigate transport coefficients of a dense sym- 
metrical Lennard- Jones mixture that were calculated along a path towards a liquid-liquid 
miscibility gap ending at a coexistence state about 15% below the critical point. The main 
result of our study is that the bulk viscosity ijb increases significantly near the coexistence 
state whereas the shear viscosity r) s does not show any change near coexistence. 7] s remains 




(13) 
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to exhibit a very weak temperature dependence also when it passes the coexistence line. 
The behavior of tjb and i] s can be qualitatively understood by theories of critical dynamics 
(see Ref. (Onuki, 2002)). 

In future studies we plan to compute the transport properties also closer to the critical 
point. Of course, in such studies much larger system sizes than those used in this work have 
to be considered. Moreover, the emergence of critical slowing down requires simulations on 
much longer time scales. 
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FIG. 1: Phase diagram of the symmetrical Lennard-Jones mixture for four choices of N as indi- 
cated. The crosses at x-q = 0.10375 mark the states for which the structure and dynamics was 
investigated (note that also T = 3.0 and T = 6.0 were studied). The solid lines are fits with Eq. © 
and the dashed lines are just guides to the eye. 
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FIG. 2: Number-number density structure factor S nn (q) for the four indicated temperatures. The 
inset shows the isothermal compressibility kt as a function of temperature, kt is estimated from 
the extrapolated value S nn (q = 0) [see Eq. (JJJ)]. 
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FIG. 3: Concentration-concentration density structure factor S cc (q) for the four indicated tem- 
peratures. The inset shows the concentration susceptibility \ as function of temperature (see 
text). 
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FIG. 4: a) "Time-dependent" shear viscosity r] s (t) for the indicated temperatures. From the 
long-time plateau we read off 7] s . b) Same as in a), but now for the bulk viscosity. 
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FIG. 5: a) Shear and bulk viscosity as a function of inverse temperature, b) Ratio i]b/i] s as a 
function of temperature. 
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FIG. 6: Mean-field type extrapolation towards the "spinodal". The solid lines are fit to the data 
sets by using the functional form given by Eq. (|12j) . 
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